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ABSTRACT 

Aims. In this paper we study the dynamics of toroidal flux tubes emerging from the solar interior, through the photosphere and into the 
corona. Many previous theoretical studies of flux emergence use a twisted cylindrical tube in the solar interior as the initial condition. 
Important insights can be gained from this model, however, it does have shortcomings. The axis of the tube never fully emerges as 
dense plasma becomes trapped in magnetic dips and restrains its ascent. Also, since the entire tube is buoyant, the main photospheric 
footpoints (sunspots) continually drift apart. These problems make it difficult to produce a convincing sunspot pair We aim to address 
these problems by considering a different initial condition, namely a toroidal flux tube. 

Methods. We perform numerical experiments and solve the 3D MHD equations. The dynamics are investigated through a range of 
initial field strengths and twists. 

Results. The experiments demonstrate that the emergence of toroidal flux tubes is highly dynamic and exhibits a rich variety of 
behaviour. In answer to the aims, however, if the initial field strength is strong enough, the axis of the tube can fully emerge. Also, the 
sunspot pair does not continually drift apart. Instead, its maximum separation is the diameter of the original toroidal tube. 

Key words. Sun: magnetic fields - Magnet ohydrodynamics (MHD) - Methods: numerical 
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It is generally accepted that solar active regions, and therefore 
sunspots, are the result of magnetic flux tubes that have risen 
from the base of the convection zone and emerged at the photo- 
sphere. The emerging flux tubes are Q-shaped, at least near pho- 
tosphere, and produce bipolar sunspot pairs. This has been the 
classic al picture of sunspot formation for some time. ICowlind 
(119461) envisaged flux tubes running as girdles around the Sun 
and suggested that loops were carried upwards by convection to 
emerge as sunspot pairs. Soon after, the idea of magnetic buoy- 
ancy was put forward as the mechanism by which flux tubes 
journey from the base of the convection zone to the photosphere 
(Parker 1955; Jensen 1955). Parker demonstrated that an isolated 
magnetic flux tube in a stratified plasma under gravity must be 
buoyant, provided that it is close to thermal equilibrium with 
its suiToundings. Although the exact mechanism by which flux 
tubes rise in the convection zone is not known, most theorists 
tend to favour magnetic buoyancy. Simulations using spherical 
geometry have now been performed which show how a twisted 
flux t ube can survive rising b uoyantly though the convection 
zone dJouve & Brunll2007HFanll2b 08). 

Due to the wide variety of phenomena that exist in 
sunspots, different modelling techniques have been developed 
to describe them. These range from equilibrium models (e.g. 
Petri e & Neukirch'1999') to studies of magnetoconvection (e.g. 
Heine mann et al. 2007). A useful collection of these methods. 



and general sunspot physics, is presented in iThomas & WeissI 
d 19921 l2008h . In recent years it has become possible to sim- 
ulate the rise of flux tubes in the upper convection zone and 
their emergence into the atmosphere. The first 3D simulation of 
the emergence of a tw i sted fl ux tube in a stratified atmosphere 
was performed by|^ (l200lh . Several similar simulations have 
been performed since and are described in a recent review by 



I ArchontisI ( |2008|) . What all of these types of flux emergence sim- 
ulations have in common is that they use the same initial con- 
dition, namely a cylindrical twisted flux tube in the solar inte- 
rior. To encourage the formation of an Q-shaped loop, a density 
deficit is introduced into the tube that is proportional to e 
where y is the horizontal coordinate along the tube length and 
A is a parameter describing the length of the buoyant section of 
the tube. When the f2-loop emerges it creates a sunspot pair at 
the photosphere. Much has been gained from this model, how- 
ever, it contains some drawbacks. The cylindrical tube continues 
to emerge and the sunspot pair continually drifts apart. This is 
not seen in observati ons, where active re gions spread to a finite 
size and then decay dLiu & Zhangll2006l) . Another drawback is 
that, in simulations with this initial condition, the original axis 
of the tube remains trapped at the base of the photosphere and 
never emerges to be almost vertical when it leaves and enters 
the sunspots. In this paper we consider a different geometry for 
initial the magnetic fiel d - a toroidal flux t ube. This model was 
originally considered bv lHood et al.l (l2009h . We shall investigate 
the behaviour of emerging toroidal tubes through a range of ini- 
tial field strengths and twists. Our approach is to consider an 
idealized stratified solar interior and atmosphere. We are inter- 
ested in the interaction of the plasma and the magnetic field, 
so we solve the 3D resistive and compressive magnetohydrody- 
namic (MHD) equations. The aim is, using the toroidal model, 
to answer the following questions: 

1. Can the axis of the original flux tube emerge into the solar 
atmosphere and, if so, how? 

2. Can the two main polarities of the active region (sunspots) 
drift to a fixed distance and then stop? 

Throughout the paper, comparisons will be made with the well- 
studied cylindrical model. The outline of the paper is as follows: 
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§2 will describe the model setup and the initial conditions used 
in the flux emergence simulations. In §3 we will examine the dy- 
namics of toroidal emergence over a range of parameters. First, a 
general description of the emergence process is explained. This 
is followed by a more detailed examination of particular cases. 
§4 will summarize the conclusions. 



uniform Cartesian grid (x,y,z) of (128, 128,256) for the (non- 
dimensionalized) region -50 < x < 50, -50 < y < 50 and 
-25 < X < 85. The boundary conditions are periodic on the side 
walls of the computational box and the top and bottom bound- 
aries are closed. A damping layer is included at the top of the 
box to reduce the reflection of waves. 



2. Model setup 

In this section we outline the numerical setup and the initial at- 
mosphere and magnetic field configurations. 



2.1. Main equations 

For our numerical experiments we use a 3D ver s ion o f 
the Lagrangian remap scheme detailed in lArber et aL I (I2001h . 
This code has been used for other flux eme rgence studies 
dArchontis & Torokll2008t lArchontis et al.ll2009 l). For each time 
step, the equations are solved in Lagrangian form and are then 
remapped back onto an Eulerian grid. All of the physics is con- 
tained in the Lagrangian step. The variables are made dimen- 
sionless against photospheric values. These values are: pressure, 
Pph = 1.4 X lO'* Pa; density, pph = 3 x 10""^ kg m"-'; temperature, 
Jph = 5.6 X 10^ K and scale height H^h = 170 km. The other 



units used in the simulations are: time, fph = 25 s; speed, Mph - 
(/7ph/pph)'^^ - 6.8 km .r' and magnetic field Bph = 1.3 x 10^ G. 
The evolution of the system is governed by the following time- 
dependent and resistive (non-dimensionalized) MHD equations: 



dt 



(pu) = 0, 



-H(u- V)uj = -Vp-H(VxB)xB-H V-T-i-pg, 

— = V X (u X B) H- T/V^B, 

dt 

p + (u ■ V)ej = -pV ■ u + 7]f + gvisc, 

V ■ B = 0, 

with specific energy density 



(2) 
(3) 
(4) 



e - 



(r-i)p 



The basic variables are the density p, the pressure p, the mag- 
netic field vector B and the velocity vector u. j is the magnitude 
of current density and g is gravity (uniform in the z-direction). y 
is the ratio of specific heats and is taken as 5/3. rj is the resistiv- 
ity which is taken to be uniform. For all the experiments we take 
= 0.001. The viscosity tensor is 



T=;u Vu + Vu' --IV-U 



2.2. Initial atmosphere 

The initial stratification of the atmosphere is similar to that 
used in previous flux emergence studies (Fan 2001; Murra vet al.l 
l2006t iMacTaggart & Hood"2009). The solar interior (z < 0) is 
taken to be marginally stable to convection since in this study 
we are focussing on the emerging field. The effects of convec- 
tion are left for future work. The photosphere/chromosphere lies 
in the region < z < 10, the transition region in 10 < z < 20 and 
the corona in z > 20. The (non-dimensionaUzed) temperature is 
specified as 



m = 



1 

7.t-10)/10 



z < 0, 
0<z< 10, 
10 < z < 20, 
z>20. 



where Tcoi- - 150 is the coronal temperature. The other state 
variables, pressure and density, are found by numerically solving 
the hydrostatic equation. This gives a numerically stable equilib- 



(1) 2.3. Initial magnetic field 



As mentioned earlier, nearly all previous studies of numerical 
flux emergence use a twisted cylindrical flux tube in the solar 
interior as the initial condition. Here we present the cylinder and 
the toroidal models. 



2.3.1. Cylinder model 

The magnetic field of a twisted cylindrical flux tube (lFanll200l]) 
is given by 



(5) B = Byir)y + Bgir)0, 
where 

B,(r) = Boe-'"^'K 
Be(r) = arBy(r). 



(6) 



(7) 
(8) 



y is the direction of the tube axis and is the azimuthal direction 
in the tube cross-section, ro is the radius of the tube and = 
+ Z^. The flux tube is uniformly twisted with a denoting the 
angle of field line rotation about the axis over a unit length of the 
tube. Bo is the initial field strength at the axis of the tube. 

The pressure inside the tube differs from the surrounding 
field-free region by pdef . Balancing the radial components of the 
Lorentz force and the plasma pressure gradient gives 



d ^ ^ ^ idBi d/^def ^ 

(rBg) + + — — = 0. 

r dr 2 dr dr 



where fj. is the viscosity and I is the identity tensor The contribu- g 
tion of viscosity to the energy equation is represented by Qvisc- ~ 
The difference in the results of the simulations with and with- 
out (small) compressive viscosity is found to be negligible. The 
code accurately resolves shocks by using a combination of artifi- 
cial viscosity and Van Leer flux limiters. In such regions, heating 
is added to the energy equation. The equations are solved on a Pdef('") - B'^e^^''~^'"°{a^r^ - 2 



(9) 



Solving this gives a pressure deficit, relative to the background 
hydrostatic pressure, of 



■ 2(?r^)IA. 



(10) 
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The density deficit can then be calculated from pdef = Pdei/T(z). 
As mentioned before, to encourage an O-loop, pdef is multiplied 
by the factor e"'"^'^" to make the middle of the tube more buoyant. 
To our knowledge, all previous flux emergence studies using the 
(buoyant) cylinder model form the Q-loop using this exponential 
profile. Some studies impose a velocity profile on the cylinder 
model to ma ke it rise and emerge, rather than make the tube 
buoyant (e.g. lMagara & Loncope.2003i) . 

2.3.2. Toroidal model 



A full derivation of this model is given in iHood et al.l(l2009h so 
we shall only outline certain steps here. If, in Cartesian coordi- 
nates, the tube axis lies in the y-direction then the equations for 
the magnetic field can be written in polar coordinates {s, (p, x): 



■ zo) with s cos 4> -y and s sin c 



■Z(.h 



where zo is the base of the computational box. Under the as- 
sumption of rotational invariance, we take the magnetic field to 
be of the form 
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Fig. 1. The height-time profiles for axes with different Bq traced 
in the y = plane. Key: Bq - 1 (soHd), Bq - 3 (dot), Bq - 5 
(dash), Bq-1 (dot-dash) and Bq - 9 (triple dot-dash). 
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where A is the flux function and is constant on magnetic field 
lines. This form automatically satisfies the solenoidal constraint, 
equation (|5]l. After some manipulation, insertion of this field into 
the magnetostatic balance equation, (V x B) x B = Vpdef , yields 
the Grad-Shafranov equation 



The resulting magnetic field for a twisted toroidal tube is 
given (in Cartesian coordinates) by 



B, 
By 



Bg(r)- 



SQ 



Bg(r) — , 

s r s 



B, = B^{r)--Bg(r) . 

s r s 



(13) 
(14) 
(15) 



d^A 1 dA d^A 



ds^ s ds dx^ 



dA 



dA 



(11) 3. Parameter Study 



where - sB^. We now define a local toroidal coordinate sys- 
tem (r, 0, (p), where 

- x^ + {s - sq)^ with s - sq - r cos 6 and x = r sin 6, 

with major axis sq. Rewriting equation ( fTTT l in these local coor- 
dinates and then taking a regular expansion in the aspect ratio 
rQl SQ produces the leading order balance equation 

Bg d 1 ^Bl dpdef 

r ar 2 dr dr 

This has the same form as equation (|9]l from the cylinder model. 
We can, therefore, choose the solutions to be the same as those 
for the straight tube. These are, as in d?) and dSjl, 



Bdi 



BQe-'^l^i, Bg = arB^ = aBore-''^'''". 



Again, Bq is the axial field strength and a is the twist. The pres- 
sure difference is, again, /?def('") = B^e''^'"l'''>(a^r^-2-2a^r^)IA. 
The temperature profile is specified and the density deficit is 
given by 



Pdef 



Ble 



,11 

{a rn- 



■ 2cKV2)/(4r(z)). 



(12) 



In the simulations the entire toroidal tube is made buoyant. We 
could find more terms in the expansion but since we do not re- 
quire an exact equilibrium (as the tube will be made buoyant) we 
only consider the leading order solution. 



Now that the basic model is established, we shall investigate the 
dynamics of toroidal flux emergence by considering a range of 
parameters. In this study we look at the effects of changing the 
initial field strength Bq and the initial twist a. The other param- 
eters will be kept constant in this paper: sq - 15, vq - 2.5 and 
zo = -25. 

3. 7. Vary/ngf Bq with fixed a: generai dynamics 

In this section we consider the effects of changing the initial field 
strength, Bq, and keep the initial twist fixed at a = 0.2. This 
value is smaller than those used in previous studies (lFanll2001t 
iMurrav et al.ll2006t |Archontis & Torok 2008 ; Hood et al.ll2009l) 
and is believed to be more applicable to the Sun. We follow the 
evolution for the cases Bq = 1,3,5,7 and 9. 

From equation ( fT2] i it can be seen that the buoyancy force 
on the tube is proportional to Bq. It is then expected that tubes 
with a stronger Bq will rise faster and further than those with 
a smaller initial value. This is confirmed in the simulations and 
the height-time profiles of the tube axes are displayed in Figure 
[1] The axes are tracked by examining the change in B^ in the 
y - plane. The field line structure is also examined to confirm 
that the the axes do indeed pass through the plane at the change 
in Bx- This method is applicable since the tubes are weakly 
twisted. By rescaling the time as f = IBq, the heights reached 
by the axes are similar in the solar interior (see Figure |2|. This 
is equivalent to measuring time on the Alfven timescale rather 
than a sound timescale. Thus the heights of the tube axes are 
not only a function of time but also of initial field strength, i.e. 
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Fig. 2. The height-time profiles in the solar interior rescaled to 
time f = fBo- Key: Bq = 1 (solid), Bo = 3 (dot), Bo = 5 (dash). 
Bo = 7 (dot-dash) and Bo - 9 (triple dot-dash). 



Fig. 4. The density deficit at the tube axis (in the y - plane) 
as a fraction of the initial unsigned density deficit against height. 
Key: Bo = 5 (dash). Bo = 7 (solid). 




Fig. 3. This diagram illustrates some of the main dynamical fea- 
tures of the emergence of a toroidal flux tube. The outline rep- 
resents the outermost field line of the emerging arcade. See the 
text for more details. 



H(t) - H{tBo) where H is the height function of a tube axis. 
This agrees w ith the behaviour of cylindrical tubes as found by 
iMurrav et al.l (i2006). The two low field strength cases. Bo = 1, 3, 
do not reach the photosphere in the time the simulation is run. 
These tubes may emerge if their evolution is tracked for a much 
longer time span. The axes of the other three cases. Bo = 5, 7, 9, 
all emerge above the base of the photosphere. These cases differ 
from the cylindrical model, where the tube axis remains trapped 
near the base of the photosphere iz - 0). There is also a distinc- 
tion between the cases themselves. The magnetic field evolution 
of the strong field cases. Bo = 7, 9, differs from that of the mod- 
erate field case. Bo = 5. Their axes rise to the corona whereas the 
axis for the Bq - 5 case stops in the middle of the photosphere. 
Before we consider why this occurs we will describe the general 
behaviour of the emergence of toroidal flux tubes. 

Figure |3] displays a diagram showing some of the main fea- 
tures of the emergence process based on the dynamics found in 
the simulations. When the flux tube reaches the base of the pho- 



tosphere (z = 0), the plasma beta (/6 - /7/(|Bp/2)) decreases 
until it is 0(1). There the ma gnetic field become s subject to a 
magnetic buoyancy instability ( Murray et a l. 2006) and expands 
rapidly into the solar atmosphere which has an exponential de- 
crease in pressure with height. As the magnetic arcade expands, 
the magnitude of By decreases with height due to flux conser- 
vation. This gradient in B,. results in a Lorentz force that drives 
horizont al shear flows along the neutral line betwee n the bipolar 
sources (lManchestej200ll: [Manchester et alj|2004l) . This can be 
understood by considering the y-component of the tension from 
the Lorentz force: 

{(B ■ V)B)„ = B ■ VBy. 

The gradient of B, is negative moving in the direction of B on 
one side of the arcade and is positive moving in the direction of 
B on the other Horizontal shearing occurs from the base of the 
photosphere to the top of the arcade. 

As the magnetic field expands into the atmosphere, plasma 
drains from the top of the arcade and follows field lines down 
to the photosphere. Due to the rapid expansion of the arcade, a 
pressure and density deficit forms at the centre of the emerging 
region just above the photosphere. As the pressure at the centre 
of the arcade is smaller than that further out at the photosphere, 
the plasma, drained from the top of the arcade, flows into this 
region (see Figure[3]l. The plasma either collects there or drains 
down the legs of the toroidal tube. 

The combination of horizontal shearing and inflow can bring 
together inclined field lines and initiate magnetic reconnection. 
This process can result in the formation of new flux ropes in the 
solar atmosphere. However, the new axis will be either above or 
below the original axis, as shown in Figure |3] 

3.2. Bo = 5,7 comparison 

Now that a basic description of the dynamics of flux emergence 
has been presented we shall examine, more closely, particular 
cases that exhibit different classes of behaviour As shown before 
(see Figure [TJ the original tube axis of the Bo = 5 case stops 
rising in the middle of the the photosphere, whereas the axis of 
the Bo = 7 case rises to the corona. We shall now compare these 
cases. 
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Fig. 5. The vertical velocity, u,, (in the y = plane) as a function 
of height. Mph is the photospheric velocity scale and //ph is the 
photospheric scale height. Key: Bq - 5 (dash), Bq-1 (solid). 

Figure |4] shows how the density deficit at the axis, as a frac- 
tion of the initial unsigned density deficit, varies with height. The 
deficit is calculated by taking density values from inside and out- 
side the tube, at the same height, in the y - Q plane. This gives a 
measure of the buoyancy of the tubes. The evolution of the two 
cases is similar. Both curves rise until the tubes become neutrally 
buoyant ?Az ~ 2.5. The Bq-1 case becomes neutrally buoyant 
at an earlier time than the Bo = 5 case because the stronger field 
case rises faster. The point of neutral buoyancy corresponds to 
f SB 110 for the Bo = 5 case and f « 52 for the Bq - 7 case. 
One important point not note is that both cases become neutrally 
buoyant in the photosphere. This is not found for the cylindrical 
model where the tubes bec ome over dense before reaching the 
photosphere (iMurrav et al.ll2006) . 

Figure |5] displays how the vertical velocity, m,, at the axis 
varies with height. The Bo - 1 case achieves more than double 
the rise velocity of the Bq = 5 case in the solar interior Both 
cases initially follow a similar profile, i.e. m- increases until a 
maximum is reached and then decreases. For the Bo = 5 case, the 
velocity decreases to approximately zero at a height z = 2.4. This 
corresponds to a time of f » 110. At this height, the Bq-1 case 
has a rise velocity of m- » 0.1, at time t ^ 52 demonstrating that 
the initial choice of Bq is crucial in determining the the evolution 
of the axis properties. However, the story is more complicated as 
the axis heights in the photosphere and above are also influenced 
by draining flows. 

As previously described, plasma drains down the emerged 
magnetic arcade and then flows into a region of reduced total 
pressure {p + |B|^/2). An example of this region from the Bo = 7 
case is displayed in Figure |6]This 'square well' profile exists 
between the heights of z ~ 1.64 and z ~ 6.8 for the Bo = 5 
case and z ~ 1.64 and z ~ 5.1 for the Bq-1 case. For Bo = 5, 
the original tube axis rises slowly (compared with the Bq-1 
case) and just reaches the bottom of the pressure deficit region 
when plasma flows into it. It is this plasma that flows on top of 
the original axis and prevents its further ascent. Figures Q and [8] 
illustrate the positions of the original tube axis, for both cases, in 
relation to the field structure of the magnetic arcade at f = 100. 

As previously described, there is horizontal shearing {y- 
direction) in the arcade as it expands. This motion combined 
with inflowing plasma can lead to magnetic reconnection. 
Figures |9] and [10] shows examples of shearing and inflows for 
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Fig. 6. The 'square well' profile for the deficit in the total pres- 
sure. This deficit exists for a finite range of heights (see text) and 
plasma draining from above flows into it. This figure is for the 
Bo = 7 case at f = 100 and {y, z) = (0, 3). 




Fig. 7. The Bq = 5 case at f = 100. The original tube axis is 
represented by a red field line. Some surrounding field lines are 
traced in grey. A magnetogram is placed at the bottom of the 
photosphere (z=0) and shows B,. 



both cases, respectively. The combination of these flows does 
indeed lead to reconnection, in both simulations, and results in 
the formation of new flux ropes. In the Bo - 5 case, a flux rope 
forms above the original axis and is able to rise to the corona. 
Figure [TT] shows the field line structure of the new rope for the 
Bo = 5 case in relation to the axis of the original tube. Similar 
behavi our has been observed in simulations using the c ylinder 
model dManchester et al.l2004l;lArchontis & Torokl2008l) . In the 
Bq-1 case, the original axis emerges to the corona and a new 
flux rope forms below it. To our knowledge this has not been 
found in previous theoretical flux emergence studies that use 
the cylinder model and do not include convection. The flux rope 
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Fig. 8. The Bo = 7 case at t - 100. The original tube axis is 
represented by a red field line. Some surrounding field lines are 
traced in grey. A magnetogram is placed at the bottom of the 
photosphere (z=0) and shows B^. 
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Fig. 9. Shearing profiles at f = 100 and (y, z) = (0, 3) for Bo = 5 
(dash) and Bo = 7 (solid). 



forms directly below the original axis and the reconnection that 
creates it produces an upflow. This upflow gives an extra kick to 
the rising of the original axis and explains the steep increase in 
the u, curve for Bq = 7 in Figure|5] The field line structure of the 
new flux rope in the Bo = 7 case is shown in Figure [T2l 

In the emergence process not all the flux is transported into 
the atmosphere, some remains in the solar interior and at the 
base of the photosphere. To quantify how much flux emerges 
and how much does not we consider the horizontal flux, through 
the central y - Q plane. 



(t>h(y = 0) 
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Fig. 10. Inflow profiles at t = 100 and iy,z) = (0, 3) for Bo = 5 
(dash) and Bo = 7 (solid). There are clear inflow profiles for 
both cases. However, when the flows meet in the centre of the 
pressure deficit zone, they produce more complex behaviour 




dxdz. 



Fig. 11. The Bo = 5 case att- 100. The red field line represents 
the original tube axis. The green field line represents the axis 
of a new flux rope. The surrounding field line structure at the 
new axis is demonstrated by some field lines traced in grey. The 
original axis is pinned down in the photosphere whereas the new 
rope is at the base of the corona. The magnetogram is at z = 0. 



This integral is calculated for the regions above and below the 
base of the photosphere (z = 0). These values are shown in 
Figures [T3] and [14] as percentages of the initial O/,, through time, 
for Bo = 5 and Bo = 7. As described above, the Bq-1 case rises 
faster and emerges before the Bo = 5 case. Also, the stronger 
field case transports more flux into the atmosphere, as expected. 
The horizontal flux remains constant, for both cases, when the 
flux tube rises in the interior and has not yet reached the pho- 
tosphere. When the tube reaches the phototsphere and becomes 
subject to the buoyancy instability, horizontal flux is transported 
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Fig. 12. The Bq - 1 case at f = 100. The red field line repre- 
sents the original tube axis. The green field line represents the 
axis of a new flux rope. The surrounding field line structure at 
the new axis is demonstrated by some field lines traced in grey. 
Reconnection occurs where the grey field lines cross. An upward 
jet from this pushes the original axis higher. The magnetogram 
is at z = 0. 





60 r 


sz. 


50 h 


o 


40 h 


o 






30 \ 


'o 






20 \ 




10 \ 




\ 



20 40 60 80 100 
t/tph 

Fig. 13. The evolution of O/^y = 0, t)l<bh{y = 0, f = 0) x 100% 
for Bq - 5 (dash) and Bo - 7 (solid) in the atmosphere. 



into the atmosphere and the same amount is depleted in the inte- 
rior. 

3.3. Plasma drainage 

We have answered the first part of question one, which was asked 
in the Introduction, by showing that the original axis of a toroidal 
tube, with sufficiently strong Bq, can emerge to the corona. 
However, the problem of why this happens in the toroidal model 
and not in the standard cyl indrical model remains to be con- 
fronted. |Muniietal] (|2006[) carry out a parameter study for the 
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Fig. 14. The evolution of <D,,(y = 0, t)l<^h{y = 0, f = 0) x 100% 
for Bq - 5 (dash) and Bq -1 (solid) in the interior. 

cylinder model. For a tube of Bo = 9 they find the maximum 
height of the original tube axis reaches z ~ 2. The main dif- 
ference between the two models lies in the geometry. The legs 
of the toroidal model rise steeply as the whole arch of the tube 
rises to emerge. The cylinder model, however, with its exponen- 
tial buoyancy profile, kinks in the centre of the tube. The buoyant 
section increases as the tube rises higher. 

As the cylinder model emerges, plasma draining from the 
emerged arcade flows down to the phot osphere and collects in 
multiple dips where the axis is trapped jArchontis et al.ll2009l) . 
In the toroidal model, flows exist in the legs of the tubes that 
correspond to draining downflows. Figures [15] and [16] show a 
cut of Mj through one of the legs for the Bq - 1 case at f = 
40, 50, respectively. At f = 40 the tube is buoyantly rising and 
the vertical velocity in the cut is positive. By f = 50, however, 
plasma begins to drain down the arcade and also down the legs 
of the tube. There is a change in sign in the velocity of the cut 
and this corresponds to a draining downflow. The geometry of 
the toroidal model allows the plasma to drain down the legs and 
not collect in dips. It is this property that allows the original axis 
of the tube to emerge. 

To determine whether or not the axis can emerge in the cylin- 
der model, we perform a simple test. Instead of using the the 
standard buoyancy profile of exp(-y^//l^), we consider 

nexp(-/M^)-(«- 1), 

where n is a positive integer This is a generalisation of the stan- 
dard profile (n = 1) and will make the central part of the tube 
buoyant and the ends of the tube over dense. The size of the 
buoyant region is controlled by n and A. The reason for choosing 
this profile is to produce a toroidal-like geometry from the cylin- 
der model. When the experiment begins, the centre will rise and 
the ends sink, giving the required shape. For our experiment we 
choose n - 6, A - 20, Bq = 7 and a - 0.4. This strong twist 
is used to help prevent the breaking up of the tube in the so- 
lar interior (Emonet & Moreno-Insertis 1998). The base of the 
computational box is lowered to zq = -50 to allow the tube to 
develop a toroidal profile. Figure[T7]depicts the shape of the tube 
at f = 86 by showing an isosurface of |B| - 0.5 and a field line 

indicating the tube axis. 

As mentioned at the beginning of this section, [Murray et al.l 
(120061) found that the axis of a tube with Bo = 9 and a - 0.4 
rises to a maximum height of z ~ 2. In our experiment the, now 
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Fig. 15. The vertical velocity profile through a cut in one of the 
legs of the toroidal tube at f = 40 and (x,y) = (0,-13). The 
vertical velocity in the cut is positive since the tube is buoyantly 
rising. 
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Fig. 16. The vertical velocity profile through a cut in one of the 
legs of the toroidal tube at f = 50 and {x,y) = (0, -13). The sign 
of the velocity in the cut has now changed since plasma drains 
down the leg, creating a downflow. 



arched, cylinder axis has risen far beyond this. By f = 86 the 
height of the axis is z ^ 10, the base of the transition region. 
This confirms that it is the geometry of the toroidal model that 
enables the efficient draining of plasma and so allows the axis to 
emerge. 

This may help to explain w hy flux em e rgence 
studies that include c onvection (^ Cheung et al.l 120071: 
iTortosa-Andreu & Moreno- Insertis 2009) find that the axis 
of the cylindrical tube emerges. If convective flows can change 
the geometry of the tube from cylindrical to toroidal, then the 
axis can emerge as in our experiment. 

3.4. Varying a with fixed Bq 

In this section we investigate the effect of varying the initial twist 
a. We will consider the evolution of tubes with a = 0.2, 0.3, 0.4. 
We will take Bq - 1 for all these cases since, as described in 
the previous section, this value results in the tube axis emerging 



Fig. 17. The geometry of the cylinder model with the buoyancy 
profile n e.x^{-y^ I A^)-{n-\) exhibits a toroidal-like shape. This 
snapshot is at f = 86 and shows an isosurface of |B| = 0.5 in blue 
and a field line in red indicating the axis. 

with the field around the axis being vertical at the centre of the 
sunspots. All tubes have an initial minor radius ro = 2.5 and 
major radius sq - 15. 



3.4.1. Rise and emergence 

The twist of the magnetic field of a flux tube results in a tension 
force acting on the tube. At the top of the initial toroidal tube the 
density deficit at the axis is given by 



Pdef - 



Blia^l 



2) 



4r(z) 



The smaller the value of a the larger the deficit. Hence the axes 
of tubes with smaller values of a are more buoyant than tubes 
with larger values. Although the whole tube is made buoyant, 
the top of the tube is more buoyant since T(z) monotonically 
decreases with height in the solar interior Figure [18] shows the 
height-time profiles for the tube axes in the solar interior As 
expected, the lower the twist, the faster the rise. 

Once the tube reaches the photosphere, it slo ws down and 
then emerges by means of a buoyancy instability (Murra v et al.l 
|2006). To investigate how the initial twist influences the emer- 
gence, we look at how the unsigned vertical flux evolves with 
time. The total unsigned vertical flux in a plane (z - Zo) is de- 
fined by 



<l>..(zo) 



\B,\ dxdy. 



Here we will consider the a - 0.2 and a - 0.4 cases. Figure [T9] 
shows the height-time profile for the top of the magnetic field, 
i.e. the top of the tube and its expansion into the atmosphere. 
Figure |20] displays the evolution of O,, for the different twist 
cases in the plane z = 10 (base of the transition region). 
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Fig. 18. The axis height-time profiles for the cases a - 0.2 
(soHd), a - 0.3 (dot) and a = 0.4 (dash) in the solar interior 
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Fig. 19. The amount of twist affects the amount of flux that 
emerges into the atmosphere. Key: a = 0.2 (solid), a = 0.4 
(dash). 
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Fig. 20. The amount of twist affects the amount of flux that 
emerges into the atmosphere. Key: a - 0.2 (solid), a - 0.4 
(dash). 



Fig. 21. The y-separation in time of the maximum and minimum 
B, at the base of the photosphere (z=0). The higher the twist, the 
slower the separation. Key: a - 0.2 (solid), a - 0.3 (dot-dash), 
a = 0.4 (dash). 

The a - 0.2 case rises to the photosphere before the a - 0.4 
case. At the photosphere, however, the smaller twist case takes 
longer to become subject to the buoyancy instability than the 
larger one. Thi s means that both c ases emerge at approximately 
the same time. [Murray et alj (12006.) found similar behaviour for 
the cylinder model. They studied the buoyancy instability by 
considering a local analysis at the photosphere. Once emerged, 
the field of the a - 0.4 case rises faster and further than the 
a - 0.2 case. This is shown in Figure [191 Figure |20l shows the 
evolution of Oi, for a - 0.4 begins before that of a = 0.2. It also 
rises at a faster rate and by f = 94 there is double the amount of 
unsigned vertical flux compared with the a - 0.2 case. When at 
the photosphere, the field for the a - 0.2 case was weaker than 
that of the a - 0.4 case. Therefore, the field emerging into the 
atmosphere is also weaker, resulting in a weaker evolution of the 
total unsigned flux. 

3.4.2. Sunspot drift 

A simulation drawback of the cylinder model, which has become 
the standard in numerical flux emergence, is that the sunspot 
pair it produces continually drifts apart until it reaches the edge 
of the computational box. Although a density deficit is intro- 
duced into the tube to form an O-loop, the entire tube is in fact 
buoyant since the exp(-3'^//l^) profile makes the ends of the tube 
weakly buoyant. This is what causes the sunspots to drift apart. 
In the toroidal model, although the entire tube is made buoyant, 
the 'feet' of the flux tube are held at a fixed distance apart in 
the solar interior Figure |2T| shows the y-separation of the max- 
imum and minimum B, the base of the photosphere (z = 0) for 
the twist cases a - 0.2,0.3,0.4. As described above, the lower 
the twist, the faster the rise to the photosphere. Once the tubes 
reach the photosphere, the fields begin to spread. The sunspot 
separation increases linearly until a peak distance is reached 
and the y-separation remains approximately constant. This peak 
distance corresponds to the major diameter of the initial tube. 
In this example, 2io = 30. The constancy of the separation 
between the two main opposite polarities (sunspots) in an ac- 
tive regio n is often used as a criterion for the region's maturity 
dLiu & Zh ang 2006). Tubes of a higher twist will spread later- 
ally more slowly than tubes of a lower twist since a higher twist 
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produces a stronger tension force. The slopes for the linear sep- 
aration phase are estimated to be 0.833 for a - 0.2, 0.815 for 
a = 0.3 and 0.8 for or = 0.4. 

The maximum separation of the sunspots in these simula- 
tions is determined by the distance of the almost vertical flux 
tube legs at the base of the numerical box. This has consequences 
for the structure of the magnetic field in the interior. The classic 
picture of flux emergence, as described in the Introduction, con- 
siders flux tubes with long wavelengths in the convection zone, 
typically generated by m = 1,2 instabilities at the tachocline, 
where m is the longitudinal wavenumber. In the local Cartesian 
approximation near the surface, this can be represented by the 
cylinder model. The toroidal model, on the other hand, has 
vertical legs. There are two possible mechanisms for forming 
toroidal tubes in the convection zone. Either the modification of 
the cylindrical tube takes the form of that described in § 13.31 
namely the (enhanced) buoyant region of the cylindrical tube is 
spatially limited and takes effect deeper in the convection zone 
(rather than near the photosphere as in previous simulations) or 
the instability, at the base of the convection zone, involves a 
higher longitudinal mode, e.g. m > 10. Both cases would pro- 
duce toroidal shaped tubes with almost vertical legs. 

4. Conclusions 

In this paper we have carried out a parametric study of the emer- 
gence of buoyant toroidal flux tubes through the solar interior 
and into the atmosphere, via 3D MHD simulations,. By varying 
two of the parameters, namely Z?o and a, we have been able to 
investigate the general behaviour of the emerging tubes. 

Keeping a constant, the variation of the initial field strength 
produces a wealth of behaviour In the solar interior the buoy- 
ancy force is proportional to B^. Tubes with stronger Bq rise 
faster and further than those with lower values. In the solar in- 
terior the tubes exhibit a self-similar behaviour By rescaling the 
time to f = tBo, the axis height-time profile s for the different Bo 
lie on top of each other (iMurrav et al.ll2006l) . This shows that the 
axis heights of the tubes are not only a function of time but also 
of Bo. 

Once emerged, the evolution of the tubes can vary strongly 
depending on the choice of the initial field strength. A value be- 
tween Bo - 5 and B^-l gives a threshold between two general 
classes of behaviour. For the Bq - 1 case, the axis rises fast 
enough to emerge into the atmosphere before plasma draining 
from the field above flows into a pressure deficit and blocks its 
ascent. The plasma instead drains below the axis and produces 
reconnection upflows that further increase the height of the axis. 
For the Z?o = 5 case, the axis does not rise fast enough to es- 
cape plasma draining on top of it and so is pinned down in the 
photosphere. 

An advantage of the the toroidal model over the standard 
cylinder model in flux emergence is that the axis of the original 
tube is able to emerge into the atmosphere and so the field at the 
centre of the sunspots is vertical. This is made possible by the 
abiUty of the the plasma to drain down the legs of the tube in the 
solar interior In the cylinder model such draining is not possi- 
ble and plasma collects in dips along the tube axis. It has been 
shown, however, that by changing the buoyancy profile to pro- 
duce a toroidal shape, the tube axis can achieve greater heights 
in the solar atmosphere. 

Keeping Bo constant, the variation in a also produces inter- 
esting behaviour. In the solar interior it is found that tubes with 
lower twists rise faster than tubes with higher twists. Once at 
the photosphere, however, low twist tubes take longer to become 



subject to the magnetic buoyancy instability. This allows higher 
twist tubes to catch up and emerge approximately at the same 
time. 

It is found that the amount of flux emerged into the atmo- 
sphere depends on the value of a. The higher the twist the more 
flux is transported upwards. In the rise of a flux tube, the stronger 
the twist, the more preserved the tube remains. The field strength 
of higher twisted tubes is stronger than lower twisted tubes when 
they emerge. This manifests itself in the amount of flux that is 
transported into the atmosphere. 

Another feature of the toroidal model, which improves upon 
the cylinder model, is that the sunspots drift to a fixed separa- 
tion and then stop. This separation is determined by the major 
diameter of the original tube. The y-separation of the maximum 
and minimum values of B~ increases linearly until the maximum 
separation is reached, the diameter of the tube. The higher the 
twist of the the tube, the slower the separation rate. 

This study has revealed many important aspects of the emer- 
gence of toroidal tubes and there are many directions for fu- 
ture work to follow, such as the inclusion of convection into 
the model to study how it affects the emergence and eventual 
break up of active regions and the interaction of multiple loops 
to model solar atmospheric events. 
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